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ABSTRACT 

Flow visualization produces data in the form of two-dimensional images. If the optical components 
of a camera system are perfect, the transformation equations between the two-dimensional image and the 
three-dimensional object space are linear and easy to solve. However, real camera lenses introduce 
nonlinear distortions that affect the accuracy of transformation unless proper corrections are applied. An 
iterative least-squares adjustment algorithm is developed to solve the nonlinear transformation equations 
incorporated with distortion corrections. Experimental applications demonstrate that a relative precision on 
the order of 40,000 is achievable without tedious laboratory calibrations of the camera. 

1. INTRODUCTION 

The photogrammetric technique was originally developed to determine topography using aerial 
photography with “metric” cameras. Due to the advent of CCD (Charge-Couple Device) digital cameras, 
non-topographic applications of photogrammetry especially in the field of close-range photogrammetry 
have become increasingly important. 

A simple method for close-range photogrammetric data reduction with non-metric cameras was 
developed by Abdel- Aziz and Karara. 1 It establishes the Direct Linear Transformation (DLT) between the 
coordinates of image points, measured with a comparator, and the corresponding object-space coordinates. 
Since the DLT is a linear approach, it neglects nonlinear distortions introduced by real camera lenses. 
Several methods exist for geometric calibration of cameras. The most accurate approach is to perform a 
laboratory calibration as described by Snow et al. 2 This is a time-consuming process and requires specially 
fabricated test objects. A more convenient method is the in-situ calibration suggested by Bopp and 
Krauss, 3 in which the effects of distortion are determined together with other parameters of transformation 
by solving a set of nonlinear equations. The drawback of the in-situ calibration is that the solution may fail 
to converge to the best answer. A third possibility is a hybrid approach proposed by Cattafesta and 
Moore, 4 which applies the results of a less rigorous laboratory calibration to the in-situ calibration to 
improve the likelihood of a good answer. 

An algorithm using an iterative least-squares adjustment scheme has been developed to solve the 
nonlinear equations of the in-situ or hybrid calibration. This algorithm requires less CPU time to achieve 
the same level of accuracy in solutions as compared to the Levenberg-Marquardt algorithm used by other 
investigators. 4 Practical tests of these algorithms are conducted with luminescent temperature- sensitive 
paint (TSP) imaging on a three-dimensional swept-wing model tested in the Supersonic Low-Disturbance 
Tunnel at NASA Langley. 



2. MATHEMATICAL FORMULATION 


In the absence of distortion, the Direct Linear Transformation (DLT) 1 between a point (X,Y,Z) in 
object space and its corresponding image space coordinates (x,y) can be expressed by the linear fractional 
equations 


x = 


LgX + L W Y+ L n Z + l 


and 


y = 


LjX + L ( Y + LZZ + L s 
L 9 X + L w Y+L n Z + 1 


( 1 ) 


Distorted Image PointD 


These equations are based on the collinearity condition that the object point, perspective center of the lens, 
and ideal image point all lie on a straight line. Eqns. (1) can be solved directly for the 11 transformation 
parameters L, L L, , if there are at least six registration marks in the image whose object-space coordinates 
are known. 

Optical distortion of imperfect 
cameras will cause the image point to fall in 
a slightly different location. Figure 1 shows 
a schematic of the image geometry with 
distortion. To reestablish collinearity, 
correction terms Ax and Ay are added to 
Eqns. (1): 


v + Ax = 


y + Ay = 


L x X + L 2 Y + L,Z + L 4 
L g X + L W Y + L n Z + 1 

l 5 x + l 6 y + l 7 z + l 8 

L 9 X + L w Y+L u Z + 1 


( 2 ) 



Figure 1. - Image geometry with distortion. 

Considering lens distortions, 5 there are three symmetric parameters (K^K 2 ,K 3 ) for radial 
distortion and two asymmetric parameters ( P l ,P 2 ) for decentering distortion. In addition, two more affinity 
parameters (A,,A 2 ) to account for non-orthogonality and differential scaling of the sensor axes are very 
applicable to CCD cameras. The cumulative influence of these distortions can be expressed as 


Ax = x'^K^r 2 + K 2 r 4 + AT 3 r 6 ) + + 2x' 2 ) + 2 P 2 x'y' 

Ay = y'(^K x r 2 + K 2 r 4 + K 2 r 6 ^j + 2 P x x'y' + P 2 (r 2 + 2 y' 2 ) + A,x' + A 2 y' 

where (x',y') is the location of the distorted image point relative to the principal point ( x () , y 0 ) : 


v = x - x n 


and 


y =y-y 0 


(3) 


(4) 


and r is the radial distance of the distorted image point from the principal point: 


r 2 = x' 2 + y' 2 


(5) 


As shown in Figure 1, the optical axis of the lens is defined by the line perpendicular to the image plane 
that passes through the perspective center (X c ,Y c ,Z c ). The intersection point of the optical axis with the 
image plane is called the principal point that is not necessarily located at the geometric center of the image, 
although it is usually in this vicinity. The distance from the perspective center to the principal point is 



defined as the principal distance c . Note that c equals the focal length of a lens when the object is imaged 
at infinity. 

The transformation between image space and object space has 9 degrees of freedom. Three of these 
(x 0 ,y 0 ,c) specify the interior orientation of the camera. The other six are associated with the exterior 
orientation of the camera: (X C ,Y C ,Z C ) related to 3 translations and (w, </>,/<■) related to 3 rotations with 
respect to the object-space coordinate axes X, Y, and Z, respectively. Therefore, the 11 transformation 
parameters L, L L u have to fulfill the two constraints as noted by Bopp and Krauss 6 


(Lj 2 + l 2 + L3 2 )-(l 5 2 + l 6 2 + l j 2 ) + ^ + L(,Lw + LlL ^ ^ 2 + U - Lm + = o 


J z + f z + J 

^9 ^ MO ^ Ml 


L L + L L + L L ( A A + AAo + AAi)(AA + AAo + AAi) _ q 


( 6 ) 


A + Ao + Ai 


The interior-orientation parameters (x 0 ,y 0 ,c) can be computed from the 11 transformation 
parameters L, L L n by the following equations 4 
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The exterior-orientation parameters (X c ,Y c ,Z c ) and ( co,([),k ) can be obtained, if desired, using the 
relations 4 
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( 8 ) 


An iterative least-squares adjustment scheme is developed to solve the 11 transformation 
parameters L, L L u and up to 7 correction parameters (K l ,K 2 ,K 3 ,P l ,P 2 ,A l ,A 1 ), plus the principal point 
(x 0 ,y 0 ) that can be either solved as additional unknowns or computed from the 11 update transformation 
parameters L, L L u by Eqns. (7) after each iteration. The algorithm solves a set of nonlinear equations of 
Eqns. (2), (6) and (7) in the rewritten forms as 




where the first pair of equations are applied to each registration mark that represents the horizontal and 
vertical distance between the measured (x,y) and the computed (x,y) using the current guesses for the 
DLT and correction parameters. This produces 2N equations for N registration marks. The remaining pairs 
of equations are defined as the L-constraints for c, and c 2 , ^-constraints for c 3 and c 4 , and C-constraints 
for c 5 and c 6 , because these equations represent the constraints to be satisfied by the transformation 
parameters L,L L n , the principal point (x n , y 0 ) , and the principal distance c, respectively. The values of 
CjL c 6 are zero when the given constraints are satisfied. In this formulation, the maximum number of 
equations is ( 2N + 6) when all constraints are applied, and the maximum number of unknowns is 20, i.e., 
( T| ,L , , Xq , y 0 , K^,K 2 P\ , P 2 , Aj , A 2 ) . 

During iteration, the values of the unknowns are updated by adding a small adjustment to each. 
The adjustments are calculated by finding a least-squares solution to 

' dfJdL, L dfJdL n dfjdx 0 dfjdy 0 dfJdK j dfJdK 2 dfJdK 3 dfJdP, 6fJdP 2 dfJdA, dfJdA , , 

dg\! dL\ L dgJdL n dgjdx 0 dgjdy 0 dg l /dK l dgJdK 2 dgJdK 3 dg i /dP l dg l /dP 1 dgJdA l dg l /dA 2 

MOM MMM M MMMMM 

dfJdP i L &f N /dL n df N /dx 0 df N /dy 0 afJdK, df N /dK 2 dfJdK, df N /dP t df N /dP 2 df N /d\ &f N /dA 2 

d gnl dL \ L dg N /dL n dg N /dx 0 dg N /dy 0 dg N / dK { dg N /dK 2 dg N /dK 3 dg N /dP t dg N /dP 2 dg N /dA x dg N /dA 2 

dcJdL, L dcJdL n 000 0 00000 

dc 2 jdL, L dc 2 /dL n 000 0 00000 

dcJdL, L dc 2 /dL n dc 3 / dx 0 dc 3 jdy Q 0 0 0 0 0 0 0 

dcJdL , L dcJdL n dcjdx 0 dcjdy 0 0 0 0 0 0 0 0 

dcJdL, L dcJdL n dcjdx n dc 5 / dy 0 0 0 0 0 0 0 0 

dcJdL { L dcJdL n dcjdx 0 dc 6 /dy 0 0 0 0 0 0 0 0 

Iteration is stopped when the maximum adjustment falls below a given tolerance. A weight matrix can be 
included in the least-squares solution to place more or less emphasis on the constraint equations or to 
account for different uncertainties in the registration-mark locations. Setting the weight of the constraints 
or the registration-mark locations to zero removes them from the solution. All registration marks and 




constraints are weighted equally in this study. Each pair of constraints can also be turned “on” or “off’ by 
control keywords at user’s choice. 

The principal-point location computed in the solution sometimes tends to large errors. Three 
techniques have been implemented to deal with this problem. The first one is to bound the principal point 
within a specified range, but let it vary freely within that range in each iteration. The second technique is to 
keep the principal point constant, using the value from the initial guess in every iteration. The principal- 
point constraints can be included in this situation as well, in which case compliance will be sought by 
adjusting the L values alone, not the principal-point location. In the third technique, the principal point is 
no longer solved as unknown variables but computed directly from the L values. The expressions for 
dfJdLj and dg i /dL j are greatly simplified since Ax and Ay no longer contribute any terms to the partial 
derivatives with respect to the principal point. To simplify the code, when the originally described method 
is implemented, partial derivatives are first computed considering x 0 and y 0 to be additional variables. 
Then the following chain-rule operation is performed: 

-dfJdL, L dfJdL n 

dgjdh L dgi/dh , 

MO M 
df N !dL x L df N / dL u 
dg N /dL l L dg N /dL u 

x 0 and y 0 are then removed from the vector of unknowns and the d/ dx () and d/dy 0 elements are removed 
from the matrix of partial derivatives. Of course, the E-constraints are not applicable under this situation. 

3. PRACTICAL TESTS AND RESULTS 

A Photometries camera (Model TEK512AF) with an AF Micro-Nikkor 60 mm f/2.8 lens was used 
in TSP luminescence imaging of a supersonic swept-wing model. This camera is a scientific-grade, 

thermoelectrically cooled, 14-bit digital CCD camera. It has a512 x 512 sensor array with 27 am square 
sensor elements. Results of laboratory calibration of this camera were given in Reference 4 and 
summarized as follows: 

x 0 (mm ) = -0.1430 
y 0 (mm ) = 0.0917 

c(mm) = 63.9537 (11) 

Aj(l I mm 2 ) = -1.392 x 1CT 5 
2 = 28000 

where 2 is the relative precision. A relative precision of 10,000 implies that 0.1 mm can be resolved in an 
image spanning 1 m along its diagonal. The higher-order symmetric, asymmetric, and affine distortion 
parameters were set to zero in this laboratory calibration after initial trials confirmed that they were 
negligible. Forcing these terms to zero reduces the risk of problems associated with over-parameterization 
in the calibration results. These laboratory calibration results will be used as a bench mark for comparing 
the in-situ and hybrid calibrations discussed in the following. 





Figure 2. - Swept-wing model with suction panel. Figure 3. - Raw image of swept-wing model 

showing 21 targets and 3 check points. 


The same CCD camera system with the same lens settings was used for the in-situ and hybrid 
calibrations with luminescent TSP imaging on a three-dimensional swept-wing model tested in the 
Supersonic Low-Disturbance Tunnel at NASA Langley. Figure 2 shows a sketch of the swept-wing model 
partially covered with a suction panel that is denoted by a polygon filled with a group of parallel lines. 
Details on the tunnel, model and TSP technique were given in Reference 7. The TSP technique is 
developed for the purpose of transition detection in order to evaluate the Laminar Flow Control technique 
with surface suction on the swept-wing model. Figure 3 shows an original “raw” image, from which the 
transition data on the lower half of the model is obtained. Note the locations of the black targets, the 
object-space centroid locations of which were measured by a Brown and Sharpe Coordinate Measuring 
Machine in terms of model coordinates. There are 21 targets used as registration marks for camera 
calibration, and 3 check points denoted as “a”, “b” and “c” are reserved for calibration evaluation. 

Figure 4 shows a “transition” image of the 
swept-wing model in Mach 3.5 flow that is from 
left to right. The long dimension of image covers 
approximately 340 mm. The onset of transition is 
demarcated in the image by a light-gray parabolic 
band. To determine the locations of the transition 
front in model coordinates, various combinations 
of the photogrammetric techniques described 
earlier in this paper were used to perform the 
mapping. Tables 1(a) and 1(b) summarize the 
results of photogrammetric mapping by the 
iterative least-squares adjustment algorithm 
described in the previous section and the 
Levenberg-Marquardt algorithm used by other 
investigators, 4 respectively. Four cases are 


Transition 



Figure 4. - Transition image of swept-wing 
model at Mach 3.5. 

included in tables: (1) the basic DLT (Eqns. 1 only), (2) DLT with the L-constraints (Eqns. 1 and 6), (3) 
the in-situ calibration including one distortion-correction parameter K l as the only one additional unknown 
and applying the L-constraints (Eqns. 2 and 6 with application of Eqn. 10), and (4) a hybrid scheme using 



the calibrated values of (L 0 ,y 0 ,c) from Eqns. 11, by holding these parameters fixed but still allows some 
flexibility in allowing the values of appropriate distortion-correction parameters to be obtained in-situ ( Aj 
only in this example; solving Eqns. 2, 6 and 7). All computations were carried on a Sun SPARC station 2 
computer. Note that both algorithms are able to obtain similar results of camera parameters and achieve 
similar accuracy that is represented either by the relative precision 2 , estimated as 340^/2 /ma x(residual), 
or by the mapping errors at the check points a, b and c. However, the least-squares adjustment algorithm 
takes much less computer time ( CPU) to reach the solutions when iteration is required. Some interesting 
points are noted. First, the A-constraints have no significant effects on the solutions. Second, reasonable 
values are obtained for c in all cases, but the accuracy associated with the principal-point location (x 0 ,y 0 ) 
is poor except for Case 4 in which (v 0 ,y 0 ,c) are held fixed at the calibrated values. Third, the accuracy 
associated with (v 0 ,y 0 ,c), computed by Eqns. 7, is improved by including Aj as shown in Case 3. 
Because the distortion corrections are proportional to the distance from (x f) , y (l ) , raised to some power, the 
errors associated with poor estimates for the principal-point location can be significant even for a camera 
system with moderate distortion. 

Further attempts to improve the accuracy associated with the interior-orientation parameters were 
tried by solving ( x 0 ,y 0 ) as additional unknowns, including only one distortion-correction parameter Aj , in 
photogrammetric mapping by the iterative least-squares adjustment scheme. Results of various 
combinations of techniques are summarized in Table 2. The A-constraints were applied to each case listed 
in Table 2 even though they do not have significant effect on the solutions. Additional techniques applied 
to the cases of Table 2 are: (1) the P-constraints, (2) the P-constraints and a range-limit for (x n , y 0 ) , (3) the 
P-constraints and C-constraints, (4) the C-constraints and a range-limit for (x 0 , _y (l ) , and (5) the P- 
constraints, C-constraints and a range-limit for (x f) , _y 0 ) . When the C-constraints were applied, the value of 
the principal distance was given as c = 62.5394 mm that was obtained from the Case 3 of Table 1(a). The 
range-limit for (x n , _y (l ) was set as ±1 pixel-dimension, i.e., 54 um with 2x2 binning, in the 
neighborhood of the image-plane geometric center. The (x f) , y 0 ) obtained from iterative solutions are 
always fallen right on the preset range-limit when it is applied. Therefore, the values of ( x 0 ,y 0 ) listed in 
Table 2 were computed by Eqns. 7 in order to be consistent with that shown in Table 1(a). Two important 
findings must be addressed. First, the results of Case 1 in Table 2 are exactly the same as that of Case 3 in 
Table 1(a). Although their approaches are different, their procedures are mathematically identical. Second, 
application of the C-constraints and/or a range-limit for (.r 0 ,y 0 ) improves the accuracy of mapping as 
denoted by the increases of the relative precision 2 . 

Table 1. Results of Photogrammetric Analysis. 

(a) Least-Squares Adjustment Algorithm. 


Parameter 

DLT 

Constrained DLT 

In Situ w / K, 

Hybrid w / K, 

x 0 (mm) 

6.566 

6.514 

3.320 

-0.143 

y 0 (mm) 

-0.736 

-0.492 

-2.279 

0.092 

c (mm) 

58.149 

58.214 

62.539 

63.954 

® (deg.) 

87.793 

88.079 

85.460 

88.305 

(j>(deg.) 

46.295 

46.356 

49.475 

52.442 

K (deg.) 

1.624 

1.403 

3.407 

1.129 

X c (mm) 

915.563 

916.053 

955.366 

967.815 

Y c ( mm) 

482.545 

482.939 

515.506 

529.822 

Z c (mm) 

2.340 

1.854 

0.767 

-3.466 

K t ( mm 2 ) 

- 

- 

3.957xl0" 5 

5.562xl0‘ 5 

I 

57,800 

57,600 

15,800 

45,300 

a error (x,y) ( /vn) 

(-31,3) 

(-30,3) 

(-32,5) 

(-31,3) 

b error (x,y) (/.im) 

(14,24) 

(14,24) 

(15,26) 

(16,25) 




Parameter 

DLT 

Constrained DLT 

In Situ vtV K, 

Hybrid vtV K, 

x„ (mm) 

6.566 

6.514 

3.957 

-0.143 

y 0 (mm) 

-0.736 

-0.492 

-2.258 

0.092 

c (mm) 

58.149 

58.214 

62.370 

63.954 

® (deg.) 

87.793 

88.079 

85.529 

88.306 

(j>(deg.) 

46.295 

46.356 

48.868 

52.441 

K (deg.) 

1.624 

1.403 

3.342 

1.128 

X c (mm) 

915.563 

916.053 

953.821 

967.811 

Y c ( mm) 

482.545 

482.939 

514.632 

529.823 

Z c (mm) 

2.340 

1.854 

0.764 

-3.454 

Kj ( mm 2 ) 

- 

- 

4.269xl0" 5 

5.605x10 s 

X 

57,800 

57,600 

12,200 

45,200 

a error (x,y) ( i m ) 

(-31,3) 

(-30,3) 

(-32,5) 

(-31,3) 

b error (x,y) ((im) 

(14,24) 

(14,24) 

(15,26) 

(16,25) 

c error (x,y) ((mi) 

(12,20) 

(12,20) 

(13,18) 

(12,18) 

CPU (sec.) 

0.5 

4.5 

150.5 

16.2 


Table 2. Results of the In-Situ Calibration with Various Improvement Techniques. 


Parameter 

1. /pc 

2. /pc, pplim 

3. /pc, /cc 

4. /cc, pplim 5. /p 

c, /cc, pplim 

x 0 (mm) 

3.320 

5.268 

6.634 

2.410 

2.657 

y 0 (mm) 

-2.279 

-1.660 

-1.552 

-1.760 

-2.878 

c (mm) 

62.539 

60.290 

62.539 

62.293 

62.527 

co (deg.) 

85.460 

86.526 

86.529 

85.758 

84.677 

<t>(deg.) 

49.475 

47.614 

46.209 

50.390 

50.120 

k (deg.) 

3.407 

2.605 

2.488 

3.237 

4.065 

X c (mm) 

955.366 

935.379 

955.187 

952.669 

956.225 

Y c ( mm ) 

515.506 

499.265 

519.916 

512.039 

514.942 

Z c (mm) 

0.767 

2.408 

-0.687 

-2.155 

1.951 

K t ( mm 2 ) 

3.957xl0~ 5 

1.872xl0- 5 

6.628x10 5 

3.340xl0 5 

3.091x10 5 

X 

15,800 

33,800 

3,900 

42,200 

32,500 

a error (x,y) ( (mi) 

(-32,5) 

(-44,7) 

(-33,6) 

(-43,4) 

(-39,9) 

b error (x,y) ((im) 

(15,26) 

(15,26) 

(15,27) 

(13,26) 

(15,28) 

c error (x,y) ( (im ) 

(13,18) 

(16,16) 

(13,17) 

(14,15) 

(17,14) 

S'otc: /pc - with E-constraints, /cc 

- with C-constraints, 

pplim - with range-limit for (x 0 , 

To)- 


Other attempts to include more distortion-correction parameters in the in-situ calibration have been tried. 
However, all of the results indicated that K ] is still the most significant parameter and the precision of 
mapping is deteriorated when more unknowns are involved. These results are not included in the paper. 

In all cases described above, the hybridthe transition data. Before mapping the model coordinates 
scheme of Case 4 of Table 1(a) has the bestto the image plane by Eqns. 2, the optical distortion was 
relative precision 2, with the most accuratecalculated by Eqns. 3 with calibrated K x and ( x 0 , y 0 ) . 
interior-orientation parameters, (v 0 ,y 0 ,c).Then, bilinear interpolation was used to determine the 
Therefore, it is chosen to perform the finalsurface values associated with the calculated pixel 
photogrammetric mapping. Figure 5 shows thelocations. Checks are made to insure that the mapped data 
mapping results. A surface grid of the modelfall within the physical boundaries of the CCD array and 
consisting of 101x101 points was used to map the ( X,Y,Z ) locations of interest are visible. 

















Figure 5. - Photogrammetric mapping of 
transition image data to CFD grid. 


o 



4. SUMMARY 


With no prior knowledge of the interior-orientation and distortion-correction parameters, the 
iterative least-squares adjustment scheme can be used to solve the nonlinear transformation equations with 
these parameters as additional unknowns. Therefore, the camera is calibrated simultaneously with solving 
the perspective transformation. This procedure is the so-called in-situ calibration. No special calibration 
rigs are needed. The computational requirements are easily within the capabilities of modem personal 
computers, and the algorithms are efficient. The major drawback is that in some cases, the iterative 
solution may fail to converge to the best solution. There exists a strong correlation between many of the 
parameters, making it very difficult to solve for the full set of distortion corrections even under the best of 
circumstances. With only one view of the registration marks, the interior-orientation parameters are very 
difficult to accurately determine, which further hinders accurate calibration. To circumvent this obstacle, a 
hybrid scheme is proposed to incorporate the laboratory-calibrated values of interior-orientation parameters 
in the in-situ calibration. A relative precision of 1 part in 45,300 is achieved in a sample application to 
luminescent temperature-sensitive paint imaging. However, when there is no laboratory calibration 
available, various techniques are proposed to improve the accuracy of the in-situ calibration. The best 
relative precision that can be achieved by the improved in-situ calibration in the same sample application is 
1 part in 42,200, that is comparable to the hybrid calibration. This suggests that the in-situ calibration 
incorporated with the iterative least- squares adjustment scheme is a viable tool for photogrammetric flow 
visualization. 
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